#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <math.h>

void seqScan(double *X, unsigned int n, unsigned int offset, 
		unsigned int stride);

void seqAdd(double *X, double val, unsigned int size, unsigned int offset, 
		unsigned int stride);

// blockBits is the number of bits for a block, presumable fits in cache.
void hybridScan(double *X, int n, int blockBits) {
	int i, j, d, logn, step, limit, mask, bit;

	logn = (int)ceil(log2(n));

	// shift the 1 into the blockBits+1 position
	step = 1 << blockBits;  
	int st;
	
	// first do blockwise sequential scan
#pragma omp parallel for shared(X, step, n) private(st) schedule(guided)
	for (i = 0; i < n; i += step) {
		st = (i + step) > n ? n - i : step;
		seqScan(X, st, i, 1);
	}
	
	// next do a parallel scan on only the "splitter" elements
	mask = step - 1;
	bit = step;
	
	for (d = blockBits + 1; d <= logn; ++d) {
		mask |= bit;  // increase mask
		step = bit << 1;

		// has to have at least mask to enter if.
#pragma omp parallel for shared(X, mask, bit, n) private(j) schedule(guided) 
		for (i = mask; i < n; i += step) 
				X[i] += X[i^bit];	// now do the math
		bit = step;
	}

	mask = (1 << logn) - 1;
	bit = 1 << (logn-1);

	// skip top bit.  each run updates 1 below.	
	// since each run udpates 1 lower in bit level,  
	// d == blockBits + 1 would update d = blockBits, 
	// which would be used to udpate the other entries in blockBits.
	for (d = logn-1; d > blockBits; --d) {  
		mask ^= bit;  // keep lower bits only
		step = bit;
		bit >>= 1;   // use to find the half step merge target.
		limit = n - bit;
		
		// update i, i is the paired target.  j has the computed sum that
		// needs to be added to i. e.g. i = 3, and j = 2.  this updates value
		// of i.

		// has to have at least mask, and can't go over n-bit.
#pragma omp parallel for shared(X, limit, n, step, bit, mask) \
		private(j) schedule(guided) 
		for (i = mask; i < limit; i += step) 
			X[i+bit] += X[i];
	}

	step = 1 << blockBits;  // shift the 1 into the blockBits+1 position	
	// finally do blockwise sequential update at the very end.
	#pragma omp parallel for shared(X, step, n) private(st) schedule(guided)
	for (i = step; i < n; i += step) {
		st = (i + step) > n ? n - i : step - 1;
		seqAdd(X, X[i - 1], st, i, 1);  // and use the previous value.
	}
}
